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Abstract. The behavior of a charge shuttle under a pure AC field has been recently 
considered theoretically and experimentally. If the system presents an asymmetry 
in the tunneling amplitudes the device acts as a nano-clectromechanical rectifier, 
transforming a pure AC voltage field into a direct currcn. In this paper we first 
review the model and the appearance of the rectifying effect for bias voltages below 
the threshold of self-oscillation. We discuss in some details the dynamics of the 
central island that, like the current, presents strong dependence on the forcing AC field 
frequency. In presence of both a constant and a small oscillating bias voltage we analyze 
the transition from the static to self-oscillating solution. We then consider current 
fluctuations (full counting statistics) for periodic motion of the grain. We explicitly 
evaluate the current noise numerically and we find that it shows clear signatures of 
correlated transport at certain locking frequencies. In the adiabatic limit we obtain a 
simple expression for the full-counting statistics and calculate explicitly the first four 
moments. 
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1. Introduction 

The field of Nano-Electro-Mechanical (NEMS) devices has been growing very rapidly 
in recent years both experimentally and theoretically Important experimental results 
in this area include the use of a Single Electron Transistor (SET) as displacement 
sensor [I] or the study of quantum transport through suspended nanotubes [3 El El, 
oscillating molecules EH IB] and islands jU HU]. On the theoretical side, several 
works [HlliaiinilIlliniliniinilIHllinillIliniE21ES! have highlighted various aspects 
of the role of Coulomb blockade in NEMS. 

An exciting prototype example of mechanical assisted SET device is the charge 
shuttle proposed in 1998 by Gorelik et al. jTJ. As shown in Ref. fl] a SET with 
an oscillating central island can shuttle electrons between the two electrodes with a 
corresponding low noise level [121 EH EI] • Although the experimental realization of the 
charge shuttle is difficult, promising systems as C 60 molecules in break junctions IE] 
or silicon structures 0111]] already exist. 

Several interesting features in the dynamics of the shuttle emerge when it is driven 
by an AC bias due to the interplay between the internal frequency of the oscillating 
island and the external frequency [21]. The most striking is probably a rectification 
effect when the contact resistances to the reservoirs are asymmetrical. Moreover since 
the whole shuttling phenomenon is due to the non-linear dependence of the electron 
transition rates on the position of the island one expects that the mechanical motion can 
respond also at multiple harmonics both of the external and of the internal resonating 
frequency. 

In the present paper we expand our results on the AC-driven shuttle presented 
in (21] and discuss in some details the behavior of the average current and its fluctuations 
together with a detailed analysis of the dynamics of the island. The paper is divided 
as follows. In Section E we present the model for the AC forced shuttle. In Section El 
we describe the dynamical behavior of the island and the current dependence on the 
external frequency. In Section 0] the transition from a static to an oscillating solution 
is studied in presence of a small AC field. In Section |S] we derive an expression for 
the FCS in the adiabatic limit where we also discuss the first four moments of current 
fluctuations. Section is devoted to the conclusions 

2. The model for the AC-driven charge shuttle 

The charge shuttle, shown in Figure ^ left panel, consists in a SET where the central 
island can oscillate between the two leads [TT] and subjected to an elastic recoil force, a 
damping force due to the dissipative medium, and an electric force due to the applied 
bias. The island is connected to the left and right leads through tunnel junctions 
characterized by resistances Rl(x) = R,L(0)e x / x and Rr(x) = Rn(0)e~ x / x . Here A is 
the tunnelling length of the order of one nanometer, and x is the displacement from the 
equilibrium position in absence of external drive. The island is coupled to the leads and 




Figure 1. Left panel: Electric scheme of a charge shuttle. Central Panel: Convention 
for the rates of charge transfer. Right panel: energy diagram for the SET. The diagonal 
lines indicate the thresholds for the vanishing of the four rates Tpp, Tpn, Ttl, and 
Ttr- The two dots indicate the state of the system during shuttling at positive voltage 
V. 



the gate through two capacitances C g and C. While it is crucial to keep into account the 
(exponential) dependence of the tunneling resistance on the position of the island, the 
dependence on x of the capacitances is much weaker and it can be ignored. The system 
is symmetrically biased at a voltage V(t) (Vr = —V L = V/2), the charging energy is 
Eq = e 2 /2C*s where e is the electron charge, Cs = C g + 2C, and the gate charge Q g is 
CgVg. We will consider the case of low temperatures {k B T <C E c ), charge degeneracy 
{Qg — e /2), and voltages \V\ < Ec/e. In this regime the grain can accommodate only 
n = or 1 additional electrons (see right panel of Figure HJ). We also assume that typical 
driving frequencies u are small compared to k R T jTi <C Eo/h. This condition on the 
frequency is rather weak (typically Ec/h ~ 10 THz) but it allows a simple description 
of the tunnelling in terms of time dependent rates. 

In the simplest approximation the dynamics of the central island can be derived 
from the Newton equation [TT| 

x{t) = -uo 2 x{t)-7 1 uj x{t) + ^^-n{t) . (1) 

together with a stochastic equation governing the charge —en(t) on the island. In 
Eq.([Tj) m is the mass of the grain, u is the oscillator eigenfrequency, t]uj is a damping 
coefficient (rj 1), and L is the distance between the two leads. In the regime of 
incoherent transport the (stochastic) evolution of the charge —en(t) is governed by the 



following four rates [25] : 

T FL = \eV(t)/4E c \T L (x)Q(V) (2) 

Tfr = \eV(t)/4E c \ T R (x) Q(-V) (3) 

Ttl = \eV(t)/4E c \ T L (x) &(-V) (4) 

T TR = \eV(t)/4E c \T R (x)Q(V) (5) 

where Eqs. (J2HH1) refers to n = 1 — ► transitions while Eqs. flUE!) to n = — ► 1 



transitions. Here EL, ER, TL, and TR stands for From/To and Left/Right, indicating 
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the direction for the electron tunneling associated to the corresponding T as shown in the 
central panel of Figure [T] T L / R (x) = [Rl/r^C] -1 and 6(t) is the Heaviside function. 
The current I is then determined by counting the net number of electrons that have 
passed through the shuttle during a given time t. 

3. Current and mechanical response to the AC field 

The problem can be studied numerically by simulating the stochastic process governed 
by the four rates Eqs.(j2HHl) an d by the deterministic evolution Eq. (JTJ). We used a 
simple technique that consists in solving analytically Newton equations between two 
tunneling events. One divides the trajectory into small intervals where the probability 
of occurrence of two tunneling events is much smaller than one, and uses the rates to 
generate events with the appropriate statistics. In the following we present results for 
typically 10 6 tunneling events for each plotted point. After a transient time the system 
reaches a stationary behavior. In Figure El the stationary DC current is shown as a 
function of the frequency of the external bias. The rich frequency landscape structure is 
generic. We observed a qualitatively identical behavior in a wide range of parameters for 
which the dissipation is sufficiently large to prevent self oscillation at constant voltage. 
The first conclusion one can draw is the existence of a rectification effect. The effect is 
of truly nanoelectromechanical origin and most probably it has been already observed 
by Scheible and Blick ^Oj- They considered a nanopilar forced by an AC field, and they 
observed a rectification effects in the MHz range. The structure of the main resonance 
strongly resembles the prediction shown in Figure El with a change of sign of the current 
at oj . The theory presented strictly does not apply to the experiment of Ref. fU] which 
was performed at room temperature, nevertheless a qualitative comparison seems to 
indicate that the main features are present even when the Coulomb blockade is not 
completely established (Eq ~ 80 K in that experiment). 

A second important point to emphasize in FigureElis the presence of a rich structure 
at low frequency. This is a consequence of the non-linearity of the problem, since the dips 
and the peaks appear at rational values of the ratio uj/uj d [cf. right panel of Figure Ej. 
In order to understand better the whole response it is useful to exploit an analytical 
description of the shuttle within the approximation scheme proposed in Ref. which 
we discuss below. 

3.1. Mean field approximation 

If the electric force is much smaller than the mechanical one, e = eV /(u^m\L) <C 1, 
one can average the force generated by the stochastic variable n{t) over many periods of 
oscillation. One should in fact realize that what matters are the hopping events. If the 
charge on the island does not vary the electric force is constant and thus no energy can be 
pumped on the system. It is thus the charge fluctuations that induce the changes in the 
dynamics. In order to have a sizable effect at e < 1 we need to take into account many 
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Figure 2. Current as a function of the frequency for e = 0.5, rj = 0.05, T/uj = 1, and 
Rr/Rl = 10. The result of the simulation of the stochastic dynamics (green points) is 
compared with the approximate I a (full line). In the small frequency region, enlarged 
in the right panel, several resonances at fractional values of w appear. We also show 
(red line) the analytical result from Eq. (|34|) in the adiabatic limit. The triangular dot 
indicates the numerical solution of the adiabatic equations. 



hopping events. This allows to neglect the fluctuations with respect to the average, even 
if the time scales of the hopping and of the oscillations are of the same order. Formally 
this consists in substituting into Eq. (JTJ) its statistical mean: (n(t)) = I] n nP n (t) where 
P n (t) is the probability that the occupation number of the island is n (in our case only 
n = and 1 are allowed values). 

The charge dynamics in the central island is then determined by the the coupled set 
of equation which include the dynamical equation for the position of the island together 
with the master equation j23] for the probability vector \P) = {P , Pi}. 

x(t) = -u 2 x(t)-r i u x(t) + ^-P l (t) (6) 
d t \P(t)) = -f(f)|P(f)> (7) 
The matrix T(t) is given by 

The probability is conserved, i.e. Pq + P± — 1. The rates depend on time through the 
extenal voltage V(t) and through the position of the island x(t). The instantaneous 
(average) current through the structure is 

i a (t)/e = p (t) r FL (t) - Pxit) r TL (t) , (9) 

where the subscript a in the current indicates that the fluctuations of the force acting 
on the shuttle, due to the discrete nature of the charge tunneling, are neglected. As 
shown in Ref. the shuttle instability, at constant bias, is controlled by the ratio e/rj. 
It is thus possible to assume that both e and r) are small, with their ratio arbitrary. 

We are now ready to discuss the dependence of the current on the external 
frequency. The most prominent structure, observed also in the experiments of Ref. [IT 
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Figure 3. Behavior of the position of the grain x(t)/X (red line, top row) compared 
to the forcing field (green line, top row). Probability P\{t) (central row) and current 
crossing the left junction (bottom row). The three sets of plots refers to the positive 
peak, to the vanishing value and to the negative peak in Figure [21 near the main 
resonance. 

is present at uj ~ uj , and it corresponds to the main mechanical resonance. The current 
changes sign across the resonance. The shape of the resonance can be explained as an 
effect of the phase relation between the driving voltage and the displacement of the grain. 
A simple way to prove this fact is to solve Eq. for a given function x(t) = A sm(ut— 0) 
therefore neglecting higher harmonics. The crucial parameter is then the phase and its 
frequency dependence. We found that if we substitute the usual relation for a damped 
harmonic oscillator: <f)(uj) = arctan^^^^o;^ — a>)] we could reproduce qualitatively the 
behavior of Fig. |2] near the resonance, and in particular the change of sign. We verified 
from the numerical simulations that <p(ui) around the main resonance agrees very well 
with the simple expression of the damped harmonic oscillator given above. 

More insight can be gained by the solution of the coupled set of Eqs. (JBJ and ((Zj) with 
periodic boundary conditions in order to obtain the stationary behavior. The results for 
three cases are shown in Figure 01 The landscape near u Q of the rectified current can 
be understood by means of a following argument. The total resistance of the structure, 
Rt{x) = Rl{x) + Rr(x), has a minimum for positive x due to the asymmetry of the 
barriers. Then, even if the transport is not given by Ohm's law, the fact that Rt(x) 
has a minimum shifted from the symmetric point allows larger currents to flow when 
the shuttle is near the right lead. For 0^0 the voltage is positive when the shuttle is 
on the right, while the opposite is true for m ir. The phase shift crosses sharply from 
to 7r when u is swept through the resonance. This explain while on the left of the 
resonance the rectified current is positive and on the right it is negative. The maximum 
is due to a combination of the fast frequency dependence of the phase and of amplitude, 
that display the usual lorenzian maximum at resonance. Between these two values the 
current has to vanish for a value of the frequency that is very near u Q . 
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From the right panel of Fig. |2l it is clear that additional structures appear for 
uj ~ u) /v with v = 2,3, . . .. Our numerical simulations indicate that, except for the 
fundamental frequency, even v are favorite with respect to odd ones. As we already 
anticipated, the motion of the shuttle and the oscillating source become synchronized 
at commensurate frequencies whose ratio is \ijv (with \x positive integer). This ratio, 
known also as the winding number, can be defined as 

w = lim 9(t)/ut , (10) 

t— >oo 

where 6(f) is the accumulated angle of rotation of the representative point (x,x) in 
the phase space. In practice 6(t) /27r gives the number of oscillations performed by the 
shuttle during the time t. It is however important to choose as a representative point 
the velocity x and the acceleration x instead of the position x and the velocity x. In 
this way the angle really gives the number of "oscillations" performed by the shuttle. 
At large value of w the structure can be quite complex as it can be seen from Figure 0] 
and the above definition succeeds in giving the right number of closed cycles in each 
case. Again these plots have been obtained in the mean field approximation. At the 
main resonances for u/oj ~ the shuttle performs v nearly harmonic oscillations 
per period of the forcing field. Since different w implies different topology in the phase 
space, it is interesting to follow the transition from one winding number to an other. In 
particular the transition from w = 2 to w = 4 is shown in Figure EJ This is performed 
without passing through the w = 3 state, that appears actually as a small plateau in 
the middle of the w = 2 plateau. 

When the system is frequency locked at a winding number w, it is possible to define 
the phase shift <ft(t) = 8(t) — wut. For an accurate analysis, the value used for w in the 
definition of (p(t) must be the best fi/u approximation to w, since even corrections of 
the order of 10 -3 to the exact fi/u results would accumulate in the calculation of <fr, and 
give an arbitrary shift of the order of tt. After a transient time for perfect locking, if 
the oscillation is perfectly harmonic, <fi should no more depend on t. Thus an additional 
important quantity to analyze is the phase shift variance (A0) 2 = (</> 2 ) — ((f)) 2 . This is 
calculated by sampling 20 points per cycle over 10 3 cycles. For a given w the smaller 
is the value of A0, the better the system locks to that external frequency. A zero 
value of the A0 indicates that the oscillation is also harmonic. A periodic, but non- 
harmonic, oscillation will induce a fluctuation of the phase in time, but the fluctuations 
are correlated and A<f> will be smaller than the completely uncorrelated result. The 
numerical results for w and A0 are shown in Fig. It shows the dependence of the 
winding number as a function of the external frequency (left panel) together with the 
analysis of A0 (right panel) calculated from the stochastic simulation, green points, and 
from the average approximation, black line. The locking at rational winding numbers is 
confirmed by the presence of plateaux of decreasing width. As expected the most stable 
plateau is at w — 1: the system locks very well at this frequency. We find that this 
holds up to very high frequency in the average approximation, where w = 1 seems the 
only possible winding number. The stochastic simulation would indicate instead that 
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Figure 4. Behavior of the shuttle in the phase space for different winding numbers 
(w is indicated in the figures). The orizontal axis is the position in units of A and the 
vertical axis is the velocity in units of Xlo . One can see for instance how the transition 
from w = 2 to w — 4 can be performed without passing through the w — 3 state. 




Figure 5. Left panel shows the calculated winding number obtained either by the 
stochastic simulation (green points) and the average approximation of Eq. (JJJ (black 
points). The inset enlarges the region near the origin. Bottom panel shows A0 obtained 
with different methods: (i) average approximation (full line), (ii) stochastic simulation 
with w from Eq. I|10(l (green points), (Hi) stochastic simulation with w set equal to 1 
(red dashed line) or (iv) with w = 2 (blue dashed line). Parameters arc the same as in 
Fig. El 
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for u/u > 1.3 the locking with w = 1, as defined from Eq. (jl(Jj) . is no more established. 
We verified that this upper limit to the w = 1 plateaux depends essentially linearly 
on the asymmetry Rr/ Rl. Quite generally plateaux size increases by increasing the 
asymmetry. However, by studying A0 for w = 1 in the whole frequency range, one 
actually obtains that a correlation is always present (i.e. A0 < n/y/3 ~ 1.81) even 
when Eq. (fTU|) gives a value of w different from one. The stochastic fluctuations thus 
unlock the shuttle globally, but not locally. Looking at the presence of local locking at 
other winding numbers we find that for instance w = 2 is clearly present for uo > 1 and 
reversely w = 1 is present around uj = 1/2 (see red dashed lines in the right panel of 
Fig. EJ). For global locking, only one phase variance is minimal. It corresponds to the 
"dominant" winding number. The presence of correlations of other winding numbers 
may indicate partial locking at these winding numbers (as in the region 0.6 ~< to <~ 0.9 
for w = 1 and 2) or the contribution of higher armonics of x(t) (as in the region to > 1). 



4. Instability in presence of the AC-drive 

In the previous section we have seen the response of a shuttle to a pure AC field for 
sub critical values of e/rj, i.e. when the shuttle instability is not present for static 
voltage. In this section we describe the transition from the stable to the unstable state 
in presence of a small forcing AC field. When the instability develops continuously from 
the static solution (soft excitation) one can study analytically the effect of an AC field 
by expanding around the static solution. We begin by solving the evolution equation in 
Fourier series: 

oo 

x(t)/\ = [A v an(vut) + B v coa{vut)] + B /2 (11) 

u=l 

The electric force for a static voltage can also be written as a Fourier series as follows: 

oo 

P n {t)n = Vv sin(z/wt) + D u cos{vujt)\ + D /2 (12) 

n u=l 

where C v = (oo/tt) sin(vut) J2 n Pn{t)n, and D v has an identical expression with the 
cosine substituting the sine. Plugging these expressions into Equation (JTJ one readily 
obtains the linear relation between the strength of the forcing harmonics (C and D) and 
the response of the system (A and B) for v > 1: 

( A„ \ = e / 1 - u 2 u 2 /uj 2 7]uu/uj \ ( C v \ . . 

{ B u ) {v 2 u 2 /ul - l) 2 + [ W u/u f \ -t,m/u 1 - v 2 u 2 /ul )\D V ) 1 ) 

The presence of lorentian resonances at uo v = uj /u is clearly visible in this form. As 
soon as a forcing exists at a certain u, the response will be stronger near the resonance 
frequency u /v. This is the reason why weak non-linearity induces structures in the 
rectified current at frequencies u j u, as shown in Figure ^ For a static voltage there is 
no reference phase, and thus we can set, for instance, B\ = 0. The coefficients C\ and 
D\ depend through the master equation (JJJ) strongly on A\ and more weakly the other 
coefficients for v > 1. Equation ()13|) may then be solved by non- vanishing values of A\ 
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only for a precise value of the frequency lu = lor. This is very near the bare oscillation 
value ui : 

u 2 R /u 2 « 1 - rjC 1 /D 1 (14) 

while the self-consistent equation for A\ becomes (near the resonance): 

r j A 1 = eD 1 (A 1 ) (15) 

The physical interpretation of this last equation is very simple, it represents the energy 
balance between the dissipated energy (left-hand side) and the pumped energy (right- 
hand side) during one oscillation jTTj. For small value of a = A\/\ it is possible to 
calculate analytically C\ and D\, it suffices to solve the master equation by expanding 
in the amplitude of oscillation a. In our case, defining T L (0) = T R (0) = Tu , at order 
a 3 we obtain: 

2f 



Ci/A = - 

' i + 4r 2 

and 



D 1 /\ 



AT 2 



1 + 



a 2 1 - 4r 2 

1 — > 16 

4 1 + AT 2 

2 1 - 12T 2 ' 

(17) 



4 i + 4r 2 

The energy balance Equation (fT3j) gives the condition e > e c = 7/(1 + 4r 2 )/r for the 
value a = to become locally unstable. The a 3 term in Equation (fTTj) determines the 
nature of the instability. If it is negative then a increases continuously from zero when 
e crosses the value e c . Note also that once the problem has been solved for the principal 
harmonics, the amplitude of the higher harmonics can be obtained using Equation fT3|) 
the solution of the first harmonics to calculate C v and D v for v > 1. 

In the presence of a small AC additional forcing field the equations (|12jl and ()13j) 
remain valid with the replacement e — > e[l + asin(u;t + 5)] and P n (£) — » P n (£)[l + 
a sin (cut + 5)], with a = Vac/Vdc ^ 1- We can again set Pi = 0, but we need to 
take into account the phase shift 5 of the forcing field with respect to the oscillation. 
By keeping only the principal harmonics we find that the equations for C\ and Di are 
modified by the appearance of the terms a cos 5/2 and a sin 5/2 on the right-hand side of 
Equation (|16J) and Equation (fT7|) . respectively. The two unknown quantities a and 5 are 
determined by Eqs. fTlj) and (|T5|) . For e<£ c one finds the usual result for a damped 
harmonic oscillator: a has a lorentian shape and the maximum value is ~ a/[2r](e — e c )] 
for 5 = 7r/2. According to this expression the amplitude diverges for e approaching 
e c , but its range of validity stops when the cubic term into Equation (J17j) becomes 
important. At that point the sin 5 term changes sign and the topological structure of 
the resonance changes drastically. In Figure IBlthe evolution of u(a) as a function of e is 
shown as obtained from a numerical solution of the system of equations ()14j) and (fT3j) 
where Pi and C\ include the linear term in a. 

Figure El shows the transition from a standard harmonic oscillator response (top- 
left panel), to a self oscillating system where periodic solutions occur only for a small 
range of the amplitude a and the frequency u (bottom-right panel). The evolution is 
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Figure 6. Evolution of the shuttle oscillation amplitude a = Ai/X as a function of 
the forcing frequency ui. The red and green points indicate stable and unstable values, 
respectively. The parameters for the plots are T = 1, r\ = 0.03, a = 0.002. This 
gives e c = 5f?. Different plots differ for the ratio e/e c that takes the values 1.03, 1.035, 
1.043, 1.06, 1.07, and 1.08, going from left to right and from top to bottom. The 
picture illustrates the evolution from a standard forced harmonic oscillator response 
to a self-oscillating shuttle. 

characterized by a change of the topology of the solution in the a-uj plane. Increasing 
the static electric field (e) the non linearity becomes more important and at some point 
(slightly larger than e c ) the lorentian shape is so deformed that an independent island 
develops (evident in the last two plots). In order to fully understand this evolution we 
need to study also the stability of the solution. This can be done by keeping a slow 
time dependence of the phase 8 and the amplitude a in deriving the equation of motion. 
Keeping only the first derivatives of a and 8 one obtains: 

d 1/2 77/4 \ / - r] a + eD 1 /X \ 

dw t\8 J ~ [ -rj/Aa 1/2 J \ a{uj 2 /uj 2 - 1) + ed/A J 1 J 

Written in this way it is clear that the solution of Equation ()14j) and Equation (JT5Jl gives 
stationary states (the right vector is zero) for the time evolution of the phase and the 
amplitude. Studying the real part of the eigenvalues of the linear expansion of Equation 
(118)) around these solutions allows the determination of their stability. Stable solutions 
are indicated in red in Figure El 

Therefore by increasing e the first part of the solutions to become unstable are those 
for uj far from the resonance. We cannot predict the behavior of the system in this case, 
we only can say that no periodic solution at the forcing frequency are possible. The 
stable regions restricts more and more around the resonance, and the shape also of the 
resonance changes drastically, as discussed above. In the end the only stable part of the 
diagram is a short segment (whose size is controlled by a) around the solution u>r, o,r 
that would be present if no AC field was applied. 
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Figure 7. Zero frequency current noise and Fano Factor as a function of the external 
frequency. The parameters are the same as in Fig.[2J). The dashed line in the left panel 
is the noise in the static SET. 



5. Full Counting Statistics 

As we have seen the asymmetry of the resistances allows to unveil the complex behavior 
of the shuttle by studying the current dependence on the forcing frequency. The 
correlated charge transfer in shuttles appears to be much more visible in the correlations 
of the current fluctuations. 

For the shuttle, at a constant bias voltage, the noise has been considered for the 
first time in ^2] neglecting the correlation among different cycles. Within a similar 
model, but taking into account the correlation between different cycles, the noise and 
the full counting statistics (FCS) of charge transfer was then calculated by one of the 
author [20]. The noise taking into account the quantum nature of the oscillator was 
obtained by Novotny et al. j2J by resorting to a numerical solution of generalized Bloch 
equations. An analytical approach was proposed at the same time in Ref. [2E]. The 
method to calculate the FCS in nanoelectromecanical systems proposed in [20] was 
extended to the quantum description of the oscillation in Ref. [2Z|. The FCS for a 
superconducting shuttle was also considered [2H1- Knowledge of the FCS gives a deeper 
insight on the dynamical evolution of the nanomechanical systems. For instance in Ref. 
|2D] the distribution of transmitted particles appeared to be strongly asymmetrical, due 
to the fact that different parameters control the transfer of one or two electrons per 
cycle. 

In the case of the AC-driven shuttle, an example of the behaviour of the noise as a 
function of the external frequancy is given in Fig. Additional interesting behaviour 
is expected in higher moments. The same method discussed in Ref. can be applied 
to study the statistics of current fluctuations also in presence of an AC field. 

Current fluctuations can be described by the probability Vt(n) of transferring n 
elementary charges during a measurement time t. After the seminal work of Levitov, 
Lee and Lesovik the theory to obtain the full counting statistics has been developed 
and applied to many different systems (for a review see for instance [20] )• It is easier to 
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calculate the generating function S that is the logarithm of the Fourier transform of Vt- 

Y,e m *V t (n, X )- (19) 



St(x) 



e St{x) 



e 

Using the definition fT9j) it is easy to show that the the n-th derivative with respect to 
i% of —S gives the n-th cumulant of transmitted particles. Explicitly the first two are 
the average particles n and the its variance (n — n) 2 . 

For the static SET the problem has been solved by Bagrets and Nazarov . They 
developed a technique to calculate in a very simple and efficient way the FCS for an 
arbitrary network of contacts in the Coulomb blockade regime (the so called "orthodox" 
regime, which is the one we are considering in the present paper). In Ref. j20| this 
technique was extended to take into account the motion of the grain, which reflects into 
a time dependence of the rates. It is possible to show that the generating function can 
be written as follows: 

(g|Texp|-^f x (t')^'}|P(0)), (20) 

here |_P(0)) is the vector introduced in Eq. (JJJ) to describe the state of the system, and 
(q\ = {1,1,..., 1}. The operator T x (t') is the rate matrix defined by Eq. (JBJ) where 
some elements have been modified. Explicitly if we want to count particles crossing a 
certain tunnel junction, we need to attach a e %x factor to the off-diagonal elements of 
(JSJ) that induces a transfer of charges from the leads to the grain, and a factor e~ %x for 
the elements involving transfer of particles in the opposite direction. In our case if we 
count particles at the tunneling contact with the left reservoir we obtain: 

\ -T FR - T FL e tx +T TR + T TL 

The problem is now reduced to solve the master equation (JZj) in presence of the counting 
field and coupled to the newton equation for x(t). When the solution is periodic the 
expression for the generetin function can be further simplified. In this case it is sufficient 
to integrate the master equation during one period: 



(21) 



A = Texp 



|-jf V x {t')dt'y (22) 

The generating function for the transfer of N electrons is determined by the iV-th power 
of A. For large N the eigenvalue Am of A with largest absolute value will dominate, 
and the initial conditions can be neglected. In this way the generating function becomes 
simply: 

-S N (x)/N = \nX M ( X ) (23) 

The problem is thus reduced in finding this eigenvalue. 

In the following we discuss the adiabatic limit of u u , t]uj , and T L (0) + r R (0). 
At leading order in u Eq. (|2T)|) becomes 

MX)/N= '\i(x,t)dt (24) 
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where \\(x, t) is the (right) eigenvalue of F x (t) at time t with minimum real part. The 
conservation of the number of particles implies that the real part is non negative. Using 
Eq. (|24jl we calculate as an explicit example the first four cumulants of the current 
fluctuation for the AC forced shuttle. For this purpose, we need the eigenvalue with 
minimal real part as a function of time for F defined by Eq. (jHJ). The variable x(t) is 
given by the solution of Newton equation neglecting the time derivatives: 

x(t) = ^§-n(t) , (25) 
where n(t) is given by the stationary solution of Eq. (JZJ): 
f+\ d f+\ F FL (t) + F FR (t) 

" (t) = P ' (t) = r rL it) + TMt) + rrUt) + TMt) ' (26) 

We consider the case e C 1, then x/X = ev(t)n(t) <C 1, where v(t) = V(t)/V a = sin(u;t). 
The exponentials in F can be expanded. It is convenient to separate the evolution into 
two parts, the range < ut < n with positive v, and the range tt < uot < 2ir, with 
negative v. The position satisfies at order e: 

^ = e V ~ [(1 + (3) Q(v) + (1-0) &(-v)] , (27) 

where (3 = [r L (0) — r^(0)]/[r L (0) + T R (0)] and we also introduce the average rate 
T = [Tl(0) + r^(0)]/2. The generating function at leading order in the adiabatic 
parameter 7 = eV F/(AE c uj) ^> 1 is S x = 7 J 2,r d(j)\(4>). A(0) is the eigenvalue with 
minimum real part matrix: f = f + ef 1. The explicit form of the matrix is: 

1+P -(1-/3) 





(28 



and 



f _ v\v\(l-/3) ( -(1 + 0) -(1-0) \ 
1_ 2 1 ^ (1 + P)j* 1-0 J 

+ 2 U( V) [-(l-0) -(1+0) ) ■ {29) 
The eigenvalue has the form: A = Ao + eAi. Ao is the eigenvalue with minimal real part 
of r . Ai can be obtained by first order perturbation theory: Ai = (/|fi|e), where |e) 
and \f) are the right and left eigenvectors of To corresponding to Ao- The calculation is 
straightforward, we give the result for the first four cumulants: 

Cl = ie ^0(l-0 2 ) (30) 



0' A )-e-0\l- 30 2 ) 



(31) 



C 2 = 1 (l-0 2 

C 3 = 76^0(1 - 2 )(1 - 120 2 + 150 4 ) (32) 
lb 

Ci = 7 (1 ~/ 2) [(1 + (3 2 - 90 A + 150 6 ) 
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e-(3 2 {l - 39/? 2 + 135/3 4 - 105/5 6 ) 



(33) 



C\ is the average number of electrons transmitted per cycle. It is clear from the 
expression that C\ and all odd cumulants vanish for symmetric structure (3 = 0. This 
can be verified in general since we find that X(x) for the first half cycle is equal to 
A(— x) f° r the second half. It can be interesting also to note that odd cumulants start 
at order e. This shows that odd cumulants probe a truly nano-mechanical effect, since 
e parameterizes the response of the island to the electric field. The explicit form of the 
first moment, the average current, is 



The corresponding value is shown in Figure 121 as a dashed line. The (small) difference 
with the numerical result comes from the e expansion, (e = 0.5 in the simulation) as 
is clear with a comparison with the black triangle obtained by solving numerically the 
adiabatic equations for the current. 

6. Conclusions 

In this work we analyzed the properties of the charge shuttle under time-dependent 
driving showing that the response of the shuttle is rather rich due to the non-linear 
dynamics of the grain. Rectification, synchronization, response to multiple of the 
external frequency appear in the electrical properties of the shuttle. All the results 
obtained here can be tested experimentally. Indeed a recent paper by Scheible and 
Blick [TU] already reports on some of the properties which we saw in the average 
current. Part of these results were presented elsewhere |2()|I21]. Here we provided a new 
description on the dynamical behavior of the system, the transition from the stationary 
to the self-oscillating solution, and the Full Counting Statistics in the presence of an 
AC-driving. 
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